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Abstract. Using the Coulomb gas method and standard methods of statistical 
c/3 . physics, we compute analytically the joint cumulative probability distribution of the 

extreme eigenvalues of the Jacobi-MANOVA ensemble of random matrices, in the limit 
of large matrices. This allows us to derive the rate functions for the large fluctuations 
to the left and the right of the expected values of the smallest and largest eigenvalues 
*^ ■ analytically. Our findings are compared with some available known exact results as 

^ . well as with numerical simulations finding good agreement. 
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1. Introduction 

Random Matrix Theory (RMT), the study of ensembles of matrices with random entries, 
has been actively explored since the pioneering works of Wigner [1] and Dyson [2], 
sixty years ago, and it still remains a very active and challenging area of research, 
with countless applications in mathematical physics, statistical mechanics and beyond. 
The three classical ensembles of RMT are the Gaussian (or Hermite), Wishart (or 
Laguerre) and Jacobi (or MANOVA) ensembles. The Gaussian ensemble appear in 
physics and are obtained as the eigenvalue distribution of a symmetric matrix with 
Gaussian entries. Both the Wishart and Jacobi ensembles are very common in statistics 
- where Wishart corresponding to the distribution of covariance matrices, and the Jacobi 
ensemble motivated by multivariate analysis of variance (MANOVA), see for example 
Muirhead [3] for real valued matrices. 

Let us first look at the Wishart ensemble which comprises of the product matrices 
W = X^X of an M x N matrix X having its elements drawn independently from 
a Gaussian distribution, P[X] oc exp[— ^Tr(X^X)]. The Dyson indices (3 = 1,2,4 
correspond to real, complex and quaternion entries of X respectively. In essence, 
Wishart matrices are random covariance matrices, constructed from random data 
matrices X (typically real entries in statistical application). By construction, such 
matrices are positive definite. Originally introduced by Wishart in 1928, these Wishart 
matrices have not only been extensively applicable in multivariate statistics, but have 
also been proved to be useful in areas of nuclear physics [3] , quantum chromo dynamics [5] 
and statistical physics [6J. Large deviation properties have also been widely investigated 
successfully in [3 [8] . 

The Jacobi/MANOVA matrix J is an N x N matrix that is a combination of two 
Wishart matrices, W(N, M x ) and W(N, M 2 ) in the form: 

7(N MM)— MjMgi) 

Ml ' M2) " W(N, Ml ) + W(N,M 2 ) ' 

where N, Mi and M 2 are the dimensions of the matrices Note that the definition of 
Jacobi/MANOVA matrix implies that Mi,M 2 > N. As will be seen below we will 
be interested in the limit of large matrices, i.e. N,M\,M 2 3> 1, such that the ratios 
c\ = N/Mi and c 2 = N/M 2 are kept constant. Therefore, we will use the notation 
<J{ci,c 2 ) from now on, with c\,c 2 < 1 (similar to the notation used for the Wishart 
matrices W(c)). The case C\ — c 2 — 1 (or Mi = M 2 = N) corresponds to square 
matrices and the case q — 1 = 0{N^ 1 ) (or Mi — N = 0(1)) is referred to as "almost 
square matrices". Last, note that by definition, the eigenvalues of J are positive and 
bounded from above by 1, namely by definition X min > and X ma x < 1- 

The deviation of the extreme eigenvalues of Jacobi ensembles from their expectation 
values plays a significant role in multivariate statistics, such as the Canonical Correlation 
Analysis (CCA). Suppose we have two data matrices, one NxMi, X = [xix 2 . . . xmJ and 
the other NxM 2 , Y = [yiy 2 . . . yM 2 }- The idea behind CCA is to seek linear combinations 
of a T x and b T y that are most highly correlated. Formally, this is done by determining 
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the largest latent root of the generalized eigenvalue problem: B — r 2 (A + B), where 
B = X T PX, A = X T P ± and P is the NxN orthogonal projection matrix Y(Y T Y)~ 1 Y T 
and P 1 - = I — P is its complement [9]. Another application of the convergence is 
in the standard generalization of the linear regression model to allow for multivariate 
responses, which is discussed in greater detail in Johnstone(2009) [ID] . Other interesting 
Jacobi applications include Quantum conductance in mesoscopic physics |Tl] and the 
multivariate F distribution [12] . 

During the years, it has been recognized that the real, complex and quaternion 
Jacobi random matrices are special cases of a larger family of random matrices now 
known as the /3- Jacobi ensemble, and which was discovered by [131 E]- I n that context, 
the real, complex and quaternion Jacobi matrices correspond to the Dyson index with 
= 1, 2 and 4 respectively, and they are special in that only for these cases there exists 
a matrix model that is not sparse. Many expressions describing statistical properties of 
these ensembles have been successfully derived, such as formulas for the moments for 
any matrix size N [15] . as well as a closed form expression for the spectral density Piv(A) 
for large N [9J [161 El EE] (essentially a generalization of the Marecenko-Pastur law for 
the Jacobi ensemble), namely 



/ \ ci + c 2 ^{x -(_){(+- xj 

MX) = xJT^x) ' (1) 

with so called natural boundaries 

(ci - l)(ci + c 2 ) + 2 c 2 ± \/cic 2 {ci + c 2 -l) 

c± ^ 5Tw (2) 

The indicator 1 x ^d takes the value 1 if x e D, and zero otherwise. As can be seen 
this spectral density has a compact support [£_, £+] suggesting that as N becomes large 
the eigenvalues are bounded in a non-trivial way. Also note that in general > 
and < 1, which reflect the absolute constraints X m in > and \ max < 1 mentioned 
above. Interestingly, for square and almost-square matrices £_ = and ( + = 1, and so 
the natural boundaries merge with the hard walls at and 1, which makes this limit 
somewhat different from the generic non-square case. 

More detailed spectral information about the Jacobi ensemble has been achieved 
in [19] where the expectation values of all the eigenvalues have been derived (again in 
the large N limit), including the expectations for the largest and smallest eigenvalues. 
Further more, the cumulative distribution functions (CDFs) for the extremal eigenvalues 
in terms of the hypergeometric function of a matrix argument have been derived in 
[201 [21] f° r an y Pi an d a simpler determinantal expression for the distribution of smallest 
eigenvalue has been derived in [18] for the classical ensembles. The determinantal 
representation lends itself to proving that the typical deviations of the largest / smallest 
eigenvalue from its expectation can be describes using the the celebrated Tracy-Widom 
(TW) distribution [22] as long as ci,c 2 < 1 (i.e., not for square or almost-square 
matrices), as has been demonstrated by Johnstone in [9]. Following this, the typical 
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fluctuations of X m ax can be expressed as (see Appendix A for more details) 
A max ~ C + + AT- 2 /3 [C l [1 ~ C+)]2 li,e Z fi + ( iV_4/3 ) > 



ClC 2 ( Cl + C2 -l)]V6' 

where Zp follows a TW distribution Fp [9j. A similar expression can be written for 



A m m as well (see Appendix A). This result is consistent with an exact inequality for the 



deviations around A max derived in [T7] , which also confirms the non-obvious scaling with 
N. 

As mentioned above, the distributions of the smallest and largest eigenvalues 
have numerous applications in multivariate statistics, but they are in practice difficult 
to calculate exactly, except for fairly small matrices. Moreover, knowledge of large 
deviations of the extreme eigenvalues is mostly unexplored. The main purpose of this 
paper is to provide exact large A" analytical results for these large fluctuations of the 
extreme eigenvalues \ max and X min from their mean values, which would be as simple 
as possible and practical. This is achieved by the use of two ingredients: the customary 
Coulomb gas mapping, which interprets the eigenvalues of a random matrix as charged 
interacting particle (see [23] and also [21] specifically for the Wishart-Laguerre case), and 
standard functional integration methods of statistical physics. This technique has been 
recently used to obtain analytically the large negative fluctuations of the maximum 
eigenvalue [7] (as well as for the minimal eigenvalue [25]) for the Wishart ensemble. 
Here we adopt this method for the Jacobi ensemble. We show that the logarithm of the 
probability of large fluctuations to the right/left of the mean value of A m j n /A max , for 
large N, can be written asymptotically as 

iog4 mm) (G) ~ -PN 2 &jT\(i - C-), Ci > C- (3) 
iogPi max) (C 2 )~-/3iV 2 ^ max) (C + -C2), C 2 <C + (4) 

where 3>±(x) are right/left rate functions (also known as the large deviation functions) 
that are derived explicitly with x being the main argument of the function. These 
expressions are valid for any large A" size of Jacobi ensembles and any value of j3 including 
(3 = 1 (real), 2 (complex) and 4 (quaternions) as special cases. 

Strictly speaking the Coulomb-gas method combined with the saddle-point 
integration can only yield the two rate functions just mentioned. However, thanks 
to an insightful physical argument by Majumdar and Vergassola developed for Wishart 
matrices [8] we will be able to derive the other two rate functions, as will be explained 
later. 

Some words about the organization of this paper. In section [2j we set up 
notations and provide some mathematical preliminaries and recall some known results 
for the Jacobi ensembles which serve to set up our notations for the rest of the 
paper. In section [31 we outline the Coulomb gas method followed by the saddle-point 
approximation. The Hilbert transformation is then used to reach a general expression 
of the Jacobi eigenvalue spectral density. In section HJ we derive the probability 
distributions explicitly for four different cases of constraints around the boundaries of 
the spectral density. Using the Majumdar- Vergassola approach [8] in section [5J we find 
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the other rate functions that could not be captured by the saddle-point method applied 
before. In section [6] the analytical predictions are compared to numerical simulations. 
Section [7] concludes the paper with a summary and discussion. Finally the technical 
derivations of integrals applied in section [3] are detailed in the appendices. 



2. Jacobi Ensembles: some preliminaries 

As explained in the Introduction, we focus on the Jacobi ensemble, with J{c\^c<i) = 
w( c ) o+vv(c2) ' wnere °i = W- = 1) 2), and W = X^ X of an (M x N) random Gaussian 
matrix X. For Mj > N (i = 1,2), the joint probability distribution function (jpdf) of 
the eigenvalues of J{c\,c<i), simultaneously derived in 1939 by Fisher, Girshick, Hsu, 
Mood and Roy (see Muirhead [3] page 112), can be written in closed form as 

N 

P N (X U . . . , X N ) oc J] X^ h - N+1) -\1 - A^^+D-i J] [A* - A/ . (5) 

i=l i<j 

Note that the jpdf in Equation fl5]) above has a symmetry under the exchange of 
(Mi <B- M2) and (A« <H- 1 — Aj). This symmetric feature of the distribution has 
important and useful consequences that will be discussed later. In short, it entails 
a certain equivalence between the distributions of the smallest and largest eigenvalues 
for the Jacobi ensemble. 

We would like to compute the probability -P/v(Cb C2) that all eigenvalues are in the 
interval £2] where < £i < (2 < 1 using the Coulomb gas method (note again that all 
eigenvalues obey < Aj < 1 by construction of the Jacobi matrix). Evidently, -P/v(Ci, C2) 
is also the cumulative probability that the minimum eigenvalue A m j n is greater than or 
equal to (1, and at the same time the maximum eigenvalue X max is less than or equal to 
C2, i-e-, 

Pjv(Ci, C2) = Prob[Ci < Ai < C2, Ci < A 2 < C2, • • • , Ci < Aiv < C2] 

= Prob[A min > Ci, X max < C2] • (6) 

In other words, P/v(Ci,C2) provides the joint probability distribution of the minimum 
and the maximum eigenvalue, which can be therefore expressed as 

Pn(Ci, C2) = ■■■ f 2 Pjv(Ai, • • • , X N )dX 1 ...dX N . (7) 

3. Analytics: The Coulomb Gas Approach 

Our starting point in the analytical treatment is the joint distribution of eigenvalues of 
the Jacobi ensemble ()5]) which can be presented as 

P JV (A 1 , . . . , X N ) = exp (~F{X)\ , (8) 
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where Z is a normalization constant, and the function F(X) is defined as 

/ 2 \ N ( 2\ - 

F(X) = - ^M 1 -N + l--j^ l °S^-[M2-N + l--j ~ At) 

TV 

-EE lo gi^-^i- ( 9 ) 

t=l j^i 

The quantity we are interested in is the joint probability of the extreme eigenvalues 
-Pjv(Ci) C2) given as a multiple integral by Eq. ((7j), which can be written as a ratio of two 
so-called (see below) partition functions 

P„(C^H^V..£^exp(-^(A)) , (10) 

with 

<2 /■& / ^ 



^(Ci,C 2 )= J rfAx - - - y dA*exp(-|F(A) 

with Z = Z(0, 1). Written in this form, -F(A) can be seen as the free energy of a two- 
dimensional (2D) Coulomb gas of charged particles (which repel each-other) constricted 
onto the real line with a logarithmic-logarithmic external potential. Z(0, 1) is simply a 
normalization constant, representing the partition function of a free or unconstrained 
Coulomb gas over A G [0, 1], while Z(Ci, £2) is the partition function of the same Coulomb 
gas but confined inside the segment AG [£1, (2}- 

The idea here is to evaluate the partition function, Z(d,(2) in the large N limit 
using the saddle-point method. 

Before proceeding further, we must first introduce the constrained charge density, 
which describes the density profile of the gas in the presence of two hard walls 
constraining the gas particles inside [Ci? C2] (from a probabilistic point of view this is 
simply the conditioned spectral density) 

1 N 

f(x) = p(A = x\Ci <x<( 2 ) = —J2 S ( x - - &)e(Ca - x) . 

i=i 

This allows us to write Z(d, ( 2 ) as a path integral over f(x) and its Lagrange multiplier 
f(x) that enforces its definition [23], [3, [25]. Z((i, £ 2 ) now takes the form 



where S[f(x)] is the action, given by 



S [f{ x )] = -J dx J dx'f(x)f(x')\og\x - x'\ - \ot x + -j^j-J J dxf (x)\ogx 



( P — 2 \ r 2 
- \ a 2 + ~^ffj dx f( x ) M 1 - x ) 

+^ dxf(x) hgf(x) + A (^J 2 dxf(x) 
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where A in the last term is the Lagrange multiplier enforcing the normalization of f(x) 

and we have defined the parameters ctj = ^ — 1, i = 1, 2 for convenience. 

We are interested in large Jacobi matrices, and consider the large N limit in this 

calculation. Thus we neglect terms of order 0{N~ X ) in the action S[f(x)}, and will 

consider the leading order terms in the action, namely 

K2 K2 K2 
S[f(x); (1,(2] = ~ / dx dx' f(x)f(x')\og\x - x'\ - a x / dxf (x)\ogx 

-a 2 ^ dxf(x) log(l -x)+A dxf(x)-lj, (11) 

which yields 

Z(Ci, C2) oc exp (-^S[f(x); Ci, Ca]) • (12) 

For large N, we can evaluate the leading contribution to the action (ITT]) using the 
saddle-point equation which reads SS[f(x)]/6f(x) = 

— 2 / dx'f*(x') log \x — x'\ — a.\ logo; — a 2 log(l — x) + A = , (13) 

where /*(#) if the value of f(x) at the saddle point. Differentiating Eq. f|T3|) once with 
respect to x, we obtain 

^ + = p r dx'i^i , (14) 



2x 2(1 — x) Jfr x — x 
with x G [Ci, C2] an d P denotes the Cauchy Principal Value. 

Equation (fT4"|) is a Tricomi integral equation whose solution involves the Hilbert 
transform. The general solution of equations of the type 

' / r(l ^ dt = ij)(x) 1 a<x<b 



71 J a t-x 
is given by Tricomi's theorem [26J: 



' T 2 y/{x - a)(b - X) \ J a t-X W 

where C is an arbitrary constant. In the particular case discussed here ( jl~4l) . we have 



/.(,) = \ - 1 p f VgEgB C _ 51 ) * + c I .(15) 

" 2 s/(x - C1XC2 - 1) \ 4 W-t) 2*. ' 



After various changes of variables, which are described in detail in |Appendix B 
il and gene 

f*(x; (1,(2) 



we 



obtain the final and general expression of the limiting spectral density: 



1 



27rv/(x-Ci)(C2-x) 



lx-e[Ci,C 2 ] • (16) 



Spectral Properties of the Jacobi Ensembles via the Coulomb Gas approach 



8 



Now we can substitute Eq. (|TB|) back into the saddle-point equation (1T51) to find the 



value of the Lagrange multiplier A which is explicitly derived in Appendix C It turns 
out that the integrals involved can be analytically solved using a change of variables, 
x = yA + d where A = ( 2 ~ Ci> an d we eventually arrive at the following expression for 
the action evaluated at the saddle point: 



S[f*(xy, £1,(2] = -{a 1 + a 2 + 2) I ailog + a 2 log 



+ ?y log + log y/(l - Ci)(l - C2) 

, VWEM± v / Ci(i-c 2 ) . C2-C1 , 17 , 

+ «ia 2 log- log — - — . (17) 

As a consequence of the symmetry highlighted previously in the joint probability 
distribution (J5J), the action ffT7|) too has symmetry properties. This can be seen clearly 
if we exchange (oi <-¥ a 2 ) and (£j •<->■ 1 — Ci), in which case the expression is mapped 
onto itself. This observation has two important applications. First, it serves as a 
reliable check for our formulae (including our spectral density in (fl6]) ). Second, it also 
allows us to predict that the expression for the probability of the minimum eigenvalue, 
P { ™ n \() with 

(0:1,0:2) is equivalent to that of the maximum eigenvalue, P^ nax \l — Q 

with (0-2, ai). 

4. The Constrained Spectral Density 

We must note that the solution in (fTEI) gives the general average density of eigenvalues 
in the limit of large for the ensemble of Jacobi matrices whose rescaled eigenvalues 
are restricted to x G [Ci ? C2] - From past experience with the Jacobi ensemble we expect 
the solution to the integral equation to yield a spectral density with a compact support, 
even when the distribution is not constrained at all (this is similar to the Marcenko- 
Pastur law in the Wishart ensemble - see [TIES]). These natural boundaries x G [£_, £+] 
which defines the unconstrained normalized eigenvalue distribution i.e. /*(x) > (and 
of course with f*(x)dx = 1), are obtained by imposing the following conditions, 

/*0r = C-;C-,C+) = o, 

/*(x = C + ;C-,C+) = 0. (18) 



This yields 



Oi(2 + «i + a 2 ) + 2 1 + Q! 2 ± a/(1 + «i)(l + a 2 )(l + «i + a 2 ) 



C± = — : 77TZ '— ^2 " > ( 19 ) 

(2 + oi + a 2 ) z 

which agrees with Eq. f fT9|) above, and hence shows that the known results can be 
recovered using the Coulomb gas approach. In essence this result is related to the 
expectation that the spectral density vanishes smoothly at the boundaries in the 
unconstrained case (and this is mathematically expressed by Eqs. (TT8~]) ). whereas this is 
not the case in the constrained cases described below. 
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Now that we have found the expression for the general spectral density, we will 
look at four cases of constraints on the spectral density. As will be described below, 
there are four cases that correspond to the relation between the imposed constraints 
(sometimes also referred to as walls) Ci> C2 an d the natural boundaries C->C+- I n order 
to simplify the description it would be useful to resort to the Coulomb gas picture, 
and to envisage the constrained spectral density as the density profile of a collection 
of charged particles that are put between impenetrable hard walls at Ci and (2- I n the 
absence of hard walls, the gas is localized inside [C-?C+]- This situation is not altered 
as long as the walls do not perturb the natural boundaries. When (1 > the wall 
penetrates the lower boundary. The wall is therefore wet by the Coulomb gas, and a 
new natural upper boundary appears, denoted heres by U(d) (obviously shifted to the 
right with respect to Similarly when (2 < C+ the upper wall is wet by the gas and a 
new lower boundary appears 1^2) • The last possibility is when both walls are effective 
and thus both are wet by the gas. 



4-1- Case 1: £1 < and C2 > C+ ~ ineffective constraints 

In this case the imposed constraints Ci, C2 he outside the natural compact support [£_, £+] 
and therefore they do not affect the eigenvalue distribution as the barriers are ineffective. 
In other words, the walls (d and £ 2 ) do not 'touch' its natural boundaries so nothing 
happens. This makes the same case as having (1 = 0, (2 = 1 or no barriers at all. 

Using the explicit values of the natural boundaries, £_ and ( + , found in Eq. ffT9l) . and 



some identities they obey (described in Appendix D ), we can write a simple expression 
for the unconstrained spectral density 

/.(*, C-, C + ) ^ l* 6[ c-,c + ] , (20) 

which is already known (see Eq. §T§ above, and Refs. El HE]), as it simply recovers 
the spectral density of the Jacobi ensemble. In essence, apart from re-deriving this 
expression using the Coulomb gas approach, we show here that in the large N limit 
conditioning on the eigenvalues being outside the natural support [£_,£+] does not 
affect anything, and in particular the spectral density remains the same. Figured] below 
shows a plot of the unconstrained spectral density when ai = 2 and «2 = 3. 



4-2. Case 2: (1 > £_ and ( 2 > C/ (Ci) _ an effective lower barrier 

In this case the barrier (1 becomes effective as it exceeds the lower boundary Also, 
(2 is ineffective as it remains > Z7 (Ci) > The density is shifted forwards to a new upper 
boundary Z7(Ci) leading to the formal density 

f( t rr/ .» 2 + a 1 +a 2 X ~ 2+j+aj y ufe) y/U(Cl) - X 

A(x;Ci,l/(Ci)) = ^ T^rW^b (21) 

obtained from the general expression (fTBl) by replacing (2 — > U(d). This density is 
plotted in figure [2] for ax = 2 and a 2 = 3. 
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Figure 2. The constrained spectral density for ot\ = 2 (i.e., c\ = 1/3) and c*2 = 3 
(i.e., C\ = 1/4), while having an effective barrier at d = 0.2 and ?7(Ci) = 0.8. 



Determining the precise value of U = C/(Ci) is achieved by requiring f*(U; (1,(2 
U) = 0. This leads to the following quartic equation for w = VU: 



w 



2Bw 3 + (B 2 + C 2 - l)w 2 + 2Bw - B 2 







where B = 71VC1 and C = 72-y/l — Ci ; an d where 7« = 



(22) 

-. This equation has four 
real solutions for all values of ctj and £1 - one negative and three positive ones. Since 
we are interested only in positive solutions, and furthermore inside the interval [0,1], 
we conclude that the right root for u>(Ci) = a/?7(Ci) is the largest solution of w, namely 

u(d) = 



B + ^+JB* 



2C 2 



z 



B(C 2 + 1) 



where z is the largest root of the cubic resolvent of equation (122]) [27] , namely 
z 3 + (2C 2 - B 2 - 2)z 2 + (1 - C 2 )(l + 2B 2 - C 2 )z - B 2 (l + C 



2\2 



0. 



z can be expressed explicitly by 

2C 2 -B 2 -2 
zo = 7; 



1 /3 

+ (r+ VQ 3 + R 2 ) +(r- VQ 3 + R 2 



1/3 



(23) 



(24) 



(25) 
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where R = ±[B 6 - 35 4 (1 - C 2 ) + 35 2 (1 + IQC 2 + C 4 ) - (1 - C 2 ) 3 ], 
andQ = -B 2 -C 2 ) 2 . 

An alternative approach to obtain the spectral density that does not require solving 
a quartic equation and therefore yields much simpler expressions, is to fix U, the 
upper support of the constrained spectral density and then work out the location of 
the corresponding lower wall Ci(U), using 

■(1 _ tf) 7l + 72V /(1 -f/) 7l 2 + [/ 72 2 -[/(!-[/) ' 2 



Ci(u) = u 



for \ [\ + i 2 - 7 2 + y/(l + 7l 2 - llf - 4 7 2 ) <U<1. 

Using the results above and making some calculations, we obtain the following 
expression for the distribution of X m in- 



P N (Kin > Ci) = Pn (Ci,U(d)) = exp 



C- < Ci < i 



(26) 



with the right rate function for the positive fluctuations of the minimum eigenvalue, 



$![™ m '(x), therefore being 



(min) 



(x) = - [S[C- + x, U{(- + x)} - S(C-, C+)] , x > , 



(27) 



and where the action S(d, C/(Ci)) = S[f*( x )] Ci> ^(Ci)] is given by: 
S((i,U) = 



(«i + a 2 + 2) «i log h a 2 log 



U 



+ ^ fog V^U + lOg y/tl-Wl-U) 



+ ai« 2 log 



V^(i-Ci) + \/Ci(i-^) 



log 



(28) 



2 ° 4 

We have also verified that E/(C-) = C+ as expected, and obtained the following 
expansion, describing fluctuations of U around ( + 



u(C- + x) = (+ + 



(i-C+)C+ 
(i-C-)C- 



-x' 



+ O (x 3 ) 



4(C + - C- 

4-3. Case 3: (2 < C+ an d Ci < ^(C2) ~ an effective upper barrier 



(29) 



Here, the barrier £ 2 crosses the point ( + from above and so this barrier becomes effective 
while (1 remains ineffective. In a similar manner as in case 2, the constrained spectral 
density shifts back to a lower bound value L(( 2 ), which also leads to the same quartic 
equation as in Eq. (j2"2"j) . only with B = an d C = 72 a/1 — £2- After making 

the same considerations as before we conclude that the right root for wfa) = \J L(( 2 ) 
corresponds to the smallest w inside [0, 1], and thus, 

2 



L(( 2 



B + v ^- JB 2 -2C 2 



B(C 2 + 1) 



(30) 
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where z is the largest root of the cubic resolvent Eq. (J23J of the quartic equation 
just as in the previous case. Finally, we can write down the full constrained spectral 
density as 



, , r / /- \ /~ \ 2 + «i + « 2 y/x-L (C 2 ) V L(Ca) X , 
/4x;L( i C2j,C2j = ^ 7^ — — lxe[L(Ca),C2 



(31) 



2tt — x x (1 — x) 

obtained from the general expression ( fT6l) by replacing Ci — >• £^2)- This density is 
plotted in figure [3] for oti = 2 and 02 = 3. 
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Figure 3. The constrained spectral density for ot\ = 2 (i.e., ci = 1/3) and a2 = 3 
(i.e., ci = 1/4), while having an effective barrier at £2 = 0.7 and L(C^) = 0.1009. 



An alternative approach to obtain the spectral density that does not require solving 
a quartic equation and therefore yields much simpler expressions, is to fix L, the lower 
support of the constrained spectral density and then work out the location of the upper 
wall (2(L), using 

r(T ,_ T [ (1 - L) 7 i + Wg - Zj 7i 2 + W 2 - L(l - L) n 2 
C2(L) " ^ [ (l-^ + ^i . 

for < L < 1 (l + 7? - 7 2 2 " V /(l + 7i 2 - 7 2 2 ) 2 - 47i 2 ) ■ 

The probability distribution of A max can then be expressed as 

Piv(A m ax<C2) = J P J v(^(C2),C2) = exp[-/3iV 2 $L max) (C + -C2)], < C2 < C+ , (32) 
where the left rate function for the negative fluctuations of the maximum eigenvalue, 

$ { r ax) (x) is 



(max) . 



1 



$: . ' (x 

and the action S(L,( 2 ) = S[f*(x)', ^((2), C2] is 



- [5[L(C+ - a;), C+ - x] - 5(C_, C+)] , x > , 



(33) 



6 (L, (, 2 J = - (en + a 2 + 2) \ a 1 log h a 2 log 
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9 2 

«i i—— a 2 



+ -f- log v^G + ~y log >/(l-^)(l-C2) 

, V(2(l-L) + y/L(l-C 2 ) C 2 ~^ 

+ aitta log ^ log — - — . (34) 

One can easily verify that L(C+) = C-> an d similar to Eq. ( 1291) for [/, we can obtain the 
following expansion describing fluctuations around 

(l-C-)C- 

L(( + -x)=C-- ^^If + (^) • (35) 

Note that by symmetry of the original joint distribution of the eigenvalues §5§ all the 
expressions in case 3, could have been obtained from those of case 2 by a\ ■<->■ a 2 , 
(f>l-( and also x 1 — x. 



4-4- Case 4'- fCi > C- an d (2 < U(Ci) ) or ((2 < C+ an ^ Ci > ^((2) ) - two effective 
barriers 

The two barriers become effective since both barriers (1 and (2 'touch' and cross their 
corresponding natural boundaries. In this case, the spectral density looks just like the 
general form in ( !T6|) . and is repeated here for the sake of completeness: 



f*{x; (1,(2) 



1 



Cti 1 



+ a 2 



27rv/(x-Ci)(C2-x) 

1 y/(l-Ci)(l-C 2 ) 
1 — x 



x 



(36) 



It is also depicted graphically in Fig. UJ 

Clearly, the action in this case will be in it's general form as expressed in (|17p and we 
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Figure 4. The constrained spectral density for a\ — 2 (i.e., ci = 1/3) and 02 = 3 
(i.e., ci = 1/4), with two effective barriers at (1 — 0.2 and (2 = 0.7 



get the partition function (1I2]) . Z(£i,£ 2 ) ~ exp[/3A^ 2 S'(Ci, C2)/2]. For the denominator, 
Z(0, 1) = Z(£_,£ + ) ~ exp[/3A^ 2 S'(C-, C+)/2], where we have used the fact that the 
solution for any (1 < and (2 > C+ ( e -g- ; when (1 = 0,(2 = 1) is the same as the 
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solution for Ci = (_ and C2 = C+ (this is exactly case 1 above). Thus, eventually the 
joint distribution of the extremal eigenvalues Pjv(Ci, C2) (T3~0T) becomes for large N 



where S(C-,C+) can be explicitly written as: 

v / ( 1 + «i)( 1 + 

«i + a 2 ) 



5(C-, C+) = - ("i + "2 + 2) a x log 



2 + «i + a 2 



y/(l + a 2 )(l + ai + a 2 )\ . «i 2 . «i 
+ «2 log + — log : 



2 + a,i + a 2 J 2 2 + ai + a 2 

a 2 2 , «2 , v / ( 1 + "iX 1 + "2) 

+ — — log h ol x ol 2 log ■ v 



2 2 + «i + a 2 - - ■ 2 + «i + a 2 

1 V^ 1 + ai)(l + a 2 )(l + ai + a 2 ) , , 

° g (2 + ax + a 2 )2 ' [6i) 

4-5. Summary of the joint distribution Pjv(A mm > Ci, A max < C2) 

It is useful to summarize the results derived above for the joint distribution -P/v(A m i n > 
Ci) A max < C2) i n one formula in order to highlight the general structure, and to facilitate 
applications: 

Pn (Amin > Cl) A max < C2) ~ 

'1 Ci < C- and C2 > C+ 

e -§NHS[ Mx )-MbU2]-S(C-,C + )} Ca < C+ and Cl < L ( C2 ) . 

e-!Jv»{S[/.(*);f 1 ,y(Ci)]-s(C-,f + )} Cl > C _ and Ca > u ( Cl ) 

^ e -|7V 2 {s[/»(x);C 1 ,C2]-5(C-,C+)} Otherwise 

where the ~ symbol stands for logarithmic equivalence for large N, or to be more 
precise, to order 0(N 2 ) in the exponential, with possible O(N) corrections. The various 
expressions for the action S are given by Eqs. (!i~7]) ; ( |28|) . f )34l) and ( )37l) . C/ (Ci) given by 
Eq. (E3D and L(C 2 ) given by Eq. (150]). 

5. The other rate functions - O(N) fluctuations on the other side 

The saddle-point approximation that we have used to evaluate Z(Ci,C2) is only able to 
capture the large fluctuations of A m ; n to the right of its mean value (A mm ) = C- (and 
similarly the large fluctuations of A max to the right of its mean value (A max ) = £+)• 
Fortunately, the authors of [8] came up with a beautiful physical argument to overcome 
this shortcoming and to estimate in their case the large deviations of the largest 
eigenvalue of the Wishart ensemble. In the following we will apply this method to 
the Jacobi ensemble, and in particular to the fluctuation to the right of A max as well as 
to the fluctuation to the left of Amir,. 
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5.1. [x): the right rate function of A max 

We will start with the derivation of the right rate function for the large positive 
fluctuations of A max , namely (x), defined using the following tail distribution 

4 max) (C) ~ e -^*+ max) «-^), C e [C+, i] • 

Note that the fluctuations in this case are of order O(N) while the fluctuations to the 
left of C+ are of order 0(N 2 ) - see Eq. (Tj0). In a nutshell, the reason for this difference 
between the right and left fluctuations is that for A max to be smaller than its mean value 
it requires ALL the N eigenvalues to be smaller than that mean value, while for A max 
to be larger than its mean value has no such mandatory effect on the other eigenvalues. 
Equivalently, using the Coulomb- gas language, pushing A max to the left compresses all 
the particles, while pulling A max to the right does not necessarily moves the others. The 
explains a difference in factor N between the energetic contributions of the two states, 
as explained below. 

The starting point of the method [8] to compute $^ ndx ) (x) is the energy expression 
(J9]). The Coulomb gas physics suggests that when the rightmost charge is moved to the 
right from ( + , i.e. A max — ( + ~ Oil), the sea of charges is a priori not subject to forces 
capable of macroscopic rearrangements. Following this physical picture, the right rate 
function is determined by the energy cost in pulling the rightmost charge in the external 
potential of the Coulomb gas and the interaction of the charge with the unperturbed sea. 
This energy cost for A max = ( 3> ( + can be estimated for large N using Equation ([9]) 

AE{() = -ai ln(C) - « 2 ln(l - () - 2 J In |C - A| p N (X) d\ , (38) 

where Ptv(A) is the average density of charges and the integral describes the Coulomb 
interaction of the rightmost charge with the sea. For large N, we use the expression for 
the spectral density (Tl6l) and the energy cost finally takes the form 

C+ 

A£(C) = -ailnC-c* 2 ln(l-C)-2 / In (C - C') /*(C'R' + C . (39) 

C- 

This energy cost expression is valid up to an additive constant C, which will be 
chosen here such that AE(£ = ( + ) = since our reference configuration is the one 
where A max = £+• We can now estimate the probability of such a configuration by 
Ftf(C)ocexp[-/9A£(C)/2]. 



Using the integrals ( 1C.1I) ( 1C.2I) detailed in Appendix C we can derive the explicit 
expression for the right rate function Q+^fe) = ^AE(( + + x) 



2 v c+; 2 v i-c 



- (2 + nj - n,)\u S Jl + JL+JJL ) _ n h, 
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where x G [0, (1 — C+)]> an d as mentioned above the additive constant was chosen to 
have $f ax) (0) = 0. Also, we use the notation A_ = C+ — C- ( a simplified expression for 



A_ is available in Appendix D ). 

We now check that the smallest large fluctuations predicted by our results match 
the largest typical fluctuations given by the Tracy-Widom distribution [9J (as long as 
Ci, c 2 < 1, i.e. for rectangular matrices). Expanding the rate functions $^ nax ' ) (a;) around 
x = gives: 

$ (max)^ x-+o 2 [(1 + oi)(l + q 2 )(l + oi + a 2 )] 1/A ^ 3/2 + 
3 C+(l — C+) 



and from Eqs. ( 133|) -(13 

$ (max) f s x^O 1 y/(l + ai)(l + a 2 )(l + «1 + « 2 ) ^3 , q, 4n 

" " ' " 24 [C + (i-C+)] 2 



This yields the following expression of -Pjv ^(C) f° r the smallest large fluctuations of 



ex P (-f x 3/2 (0) , Ce[C + ,i] 

exp(-||x(C)| 3 ), CG[0,C+] 



with 



p(max)/>x 

L ex P {-M\XKk)\-) , st 
[(1 + ai)(l + a 2 )(l + ai + a,)] 1 / 6 (C - C4 



x(C) = 



[ C+ (l-C + )] 2/3 ' 

where the symbol ~ in the last expression stand for log-equivalence for large N. 
This result agrees with Tracy-Widom distribution for large |%| (9j [TQj, [22] - see also 



Appendix A 



Last, for larger atypical fluctuations, we have the following asymptotic behaviour of the 
rate function 



$V ~ -yln(l-C+-a;) + 



a 2 In a 2 - 1 ln(l + a 2 ) 



1 + ai , , 1 + ai + a 2 ^ 
+ — - — m(l + ai) m(l + on + a 2 J 



+ 0(l-C+-x). 



5.#. &™ m \x): £/ie Ze/t rate function of X mm 

We now derive the analogous formula for the large fluctuations of the smallest eigenvalue 
A min , namely $_ m ^(x), which is defined using the following tail distribution 

4 min) (C) ~ e-^^K-o, Ce[o,C-]. 

Using the same ideas as above, similar steps and integrals (more specifically ( 1C.3I) ( 1C.4I) 



from Appendix C) we arrive at the following final expression for the rate function 



*«"»• W = !|ln(l-£) + f,n(l + T ^ 
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x 



A_ 



x 



(2 + a 1 + a 2 )\n[Jl + — + J—)-a 1 \n[Jl + — -\I^J — 



A 



x 



'C4 



X 



A_ 



x 



a 2 ln[ + 



Expanding the rate functions $^ nm - ) (x) around x = gives: 

$L min) (x 
and from Eqs. (|27l)-(l28l) 



1-C-V A 



(4 



(min) //rA x^p 2 [(1 + q 1 )(l + q 2 )(l + ox + a 2 )] 1/4 ^ 3/2 + Q , 5 / 2 <> 



C-(i-C- 



(min) / \ x— >0 

1 v /(l + ai)(l + a 2 )(l + 
24 [C-(l-C-)f 



■x 3 + 0(x 4 ) 



which yield the following expression of P^^iC) for the smallest large fluctuations of 

Amin • 



4 min) (0 



exp(-fx 3/2 (0), CG[0,C-] 

ex P (-^|x(C)l 3 ), Ce[C-,i] 



with 



x(C) = 



[(1 + «!)(! + a a )(l + ax + as)] 1 / 6 (C- - C) 



_£_)] 2 / 3 N 1 / 3 ' 

This result again agrees with the Tracy- Widom distribution for large \x\ QUI 



see also Appendix A 



Last, for larger atypical fluctuations, we have the following asymptotic behaviour of the 
rate function 



$L min) (x) x ~ C " - ^ ln(C- -x) + 



_ 1 + ai , _ , 
«i In «i m(l + acij 



l + a 2 , l + ai + a 2 
+ — - — m(l + a 2 ) m(l + a x + a 2 ) 



+ 0(C--x). 



6. Numerical Results 



Formulas f|T6|) and ffT9l) have been numerically checked using MATLAB and the 
agreement with the analytical results is already excellent. In Fig. 5, we plot the 
histograms of normalized eigenvalues A for different values of N, up to N = 100, with 
an initial sample of 10 6 Jacobi matrices ((3 = l,ci = 1/3, c 2 = 1/4). On top of each 
histogram, we plot the corresponding theoretical distribution given by Eq. (120]) . An 
important reason for presenting this figure is to give an idea about the convergence rate 
of the spectral distribution to is large- N limit. 

Next, we performed numerical simulations to confirm the correctness of the formulas 
for the constrained spectral densities. For this purpose we have used the so called 
tridiagonal models for /3- Jacobi ensembles [131 EH IZL|, described in the following. The 
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(a) (AT, Mi, M a ) = (5,15,20) 




(b) (N,M 1 ,M 2 ) = (10,30,40) 





(c) (N,M!,M 2 ) = (30,90,120) 



(d) (N,M U M 2 ) = (100,300,400) 



Figure 5. The red circles in the subfigures represent the normalized histograms of 
eigenvalues of Jacobi ensembles that have been left to run for 10 6 repetitions. From 
top left to bottom right, the figure shows how the theoretical limiting spectral density 
(blue dashed line) converges gradually with the corresponding simulation results (red 
circles) as we increase our N values, with a.\ = 2 (ci = 1/3) and a 2 — 3 (c% — 1/4). 
At large N = 100 in subfigure|(d)| the graphs can be seen to align precisely. 



underlying observation (proven in [13]) is that the tridiagonal positive semi definite 
matrix J = ZZ T is /3- Jacobi distributed (i.e., has exactly the same joint probability 
distribution of eigenvalues as where Z is the following upper bidiagonal N x N 
matrix 

/ c N -satcV-i 

CN-1 s> N-1 — s N-1C'n-2 
CN-2s'n-2 



\ 



\ 



-s 2 c'i 



(42) 



where all q's and c^'s are independent random variables distributed as follows: 



d ~ J Beta (f (A/a! + i), f (Na 2 + i)) 



1< i < N . 



(43) 
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jBeta(§i,§(N ai + Na 2 + l + i)), s' i = sjl-c' 2 l 1 < i < N - 1 , (44) 
where Beta(s,t) stands for the beta distribution, whose density is proportional to 



x 



\t-i 



on [0,1]. 



Figure |6] shows a comparison of various spectral densities simulated using the tridiagonal 
trick and our theoretical predictions based on the Coulomb gas picture. We used here 
complex (/3 = 2) Jacobi matrices (N = 30) with a,\ = 0.6, a 2 = 0.8. Fig. 6(a) show the 



unconstrained spectral density with a good agreement with the numerical simulations 
(we used here 10 5 samples of matrices). The other three insets describe the case with 



an effective lower barrier (\ (Fig. 6(b)), an effective upper barrier ( 2 (Fig. 6(c)), and 



two effective barriers (Fig. 6(d)). 

Next, we compare the Coulomb gas predictions for the probability function P/v(A max < 




(a) Unconstrained 




(b) Ci = 0.07 , C2 = C/(0.07) ~ 0.9431 

4r 




Figure 6. (Color online) j3 = 2,N = 30, a\ = 0.6, ct2 = 0.8 with 10 5 samples using the 
tridiagonal trick. Every Histogram represent the results of the numerical simulation, 
and the (red) thick line comes from the Coulomb gas prediction given by Eqs. ([20]) , 
(|3Il, (J2U and (Till). 



x) and the exact expression of Dumitriu [20J, using matrix hypergeometric functions. 
For completeness and convenience we quote below the exact result using the notation 
used in our paper 

P(x <T] Iffi (f N(2 + a 1 + q 2 )) if (f N + 1 - f ) fo (1+tt0 
( ] Tf (§N(2 + a x ) + 1 - I) Tf (fiV(l + a 2 )) 

X 2 F[ P) (fiV(l + ax), -§{Na 2 + 1); f JV(2 + ai) + 1 - |; xI N ) , (45) 
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where T^\c) is the multivariate gamma function of parameter > [20] . 

(oi, a2',bi, X) is the hypergeometric function of matrix argument X and parameter 
f3 > [31 [20] , and 1^ is the N- dimensional unit matrix. Ref. [21] gives an alternative 
expression for the density of the largest eigenvalue, as well as for the smallest eigenvalue. 
This expression is in practice hard to evaluate, and is not feasible for very large values 
of N. In the following, when considering numerical results, we have used the available 
MATLAB package by Koev and co-workers [28l [29] to numerically evaluate this matrix 
hypergeometric function for N = 5 and (3 = 2. Figure 7 shows a comparison of the 
exact result with our predictions. As can be seen the agreement is very good, which 
implies that N = 5 is already quite large for the large- N limit to hold in the Coulomb 
gas results. A few comments are in place here - as mentioned above, larger values 
of N are already challenging computationally for the matrix hypergeometric function. 
Furthermore, obtaining precise estimates of the probability density near the hard wall 
at £ = 1 is very challenging using the exact formula and the MATLAB package, which 
limits further the applicability of the exact formula. To test larger values of N and other 
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Figure 7. (Color online) Results for the largest eigenvalue distribution — In pj^ 8 ^ (£) 
vs C- Here N = 5, Mi = 8 and M 2 = 9 (ct\ = 3/5, a-i = 4/5), Jacobi matrices are 
complex (/3 = 2) . The large deviation functions (dashed lines) compare very well with 
the exact result of [201 HI] (red dots) . 



values of (3, we have performed simulations based on biased sampling j30j [31]. We adapt 
here the nice explanation in [31] to our needs. The Coulomb gas is simulated using the 
Metropolis algorithm. This works as follows: One proposes a small move A — > A', which 
is then accepted with probability 

p = min = min ( e -f(nVKW) ; 

or rejected with probability 1 — p. 

Now we want to study the statistics of the maximum eigenvalue given by Pj^ (C)- 
This distribution would be peaked around its average, so events in the tails, the ones 
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we are interested in, will be extremely rare. Thus, the standard Metropolis algorithm 
will provide a very poor sampling of the tails, as explained in |31j . 
To force the algorithm to explore regions of rare events, we use a trick not dissimilar to 
its mathematical counterpart: we put a barrier A* and write P$ nax ''(C) as 

The probability P^^CIC < A*) means the probability that the maximum eigenvalue 
takes the value ( conditioned on £ being smaller or equal to the barrier A*. This 
probability can be estimated readily using the Metropolis algorithm. Indeed, we simply 
start with a configuration such that A max < A* which is then updated according to the 
Metropolis algorithm, rejecting any move with A max > A*. This allows us to construct 
histograms of P^ (C|C — A*) for various values of the barrier A*. Note that for each 
value of the barrier A* the algorithm only allows us to explore a small region around (, 
that is C G [A* — 5, A*] for some small 5. Finally, since we can write: 

$r , (c) = _ logPr'(CIC<A>) + ^ ^ _ logP( C < A*) 



(3N 2 (3N 2 
where the constant K\* depends only on A* and ( G [A* - 5, A*] , we can estimate the 
derivative of the rate function $^ nax ' ) at ( only from the knowledge of Pjv (C|C ^ ^*)- 
In other words, full knowledge of K\* is not essential to estimate $ /(max )(£). 
We have implemented this procedure for /3 = 1, N = 50, Mi = 80 and M2 = 90 (the 
same ctj's as above) and the estimates for the the derivative of rate function $'( max )(£) is 
depicted in Figure E\ together with the predictions based on the Coulomb gas approach 
derived above. As can be seen, for N = 50 the agreement between the prediction and 
simulation is very good, and essentially indistinguishable visually. 



7. Summary and Conclusions 

In this work we study the large deviation functions of the extreme eigenvalue of the 
Jacobi ensemble of random matrices using the Coulomb gas approach. We obtain explicit 
expressions for the spectral densities, for the joint probability distribution of the extreme 
eigenvalues and both the right and left rate functions that describe large fluctuations 
around the mean value of the extreme eigenvalues. These expressions are expected 
to hold for large matrices, and a comparison to numerical simulations confirms this. 
Furthermore, the numerical results as well as the available exact expression suggests 
than even for fairly small values of N (such as N = 5) our results offer a reasonable 
approximation. This implies a potential important application of this work, namely 
providing a readily computable approximation for the probabilities of the extreme 
eigenvalues - expressions which may be useful in multivariate statistics. Among the 
numerous techniques require information about the extreme eigenvalues of the Jacobi 
ensemble are Canonical Correlation Analysis (CCA), the largest root test for a multiple 
response linear model, and more - as discussed earlier in the introduction and in many 
references in pi QUI ETJ . 
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Figure 8. (Color online) Results for the derivative of the large deviation functions 
of the largest eigenvalue distribution <£>'( max ) based on simulations. Here N = 50, 
Mi = 80 and Mi = 90 (ai = 3/5, a-i — 4/5), Jacobi matrices are real (/? = 1). The 
derivative of Coulomb gas prediction for the large deviation functions (dotted and 
dashed lines) compare very well with the simulations (red circles). The inset allows a 
better resolution in the restricted range. 



An interesting extension of this work could be the calculation of sub-leading 
corrections to the various distributions and rate functions calculated in this paper. 
This is important especially since these corrections round off the sharp boundaries at 
C± and describe the appearance of tails for any finite N. This may be achieved by 
using a recent method developed for Gaussian matrices, and presented in Refs. [321 133] 
where the the corresponding left/right rate functions are calculated. It is hoped that 
this method could be applied to the Jacobi ensemble. 
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Appendix A. Typical fluctuation of the extreme eigenvalues in terms of the 
Tracy— Widom distribution 

It is known [91 [10] that typical deviations of A max are related to the Tracy-Widom 
distribution [22]. In particular, the logit transform of A max , namely 

W = logit(A max ) = log (-^- 

approaches the Tracy-Widom distribution in the large N limit, after being properly 
centered and scaled, as follows 

W(N, oji, a 2 ) - /J>w(N, ai, a 2 ) n^oo „ 
&w{N, ai, a 2 ) 

where Fg is the TW distribution with the corresponding Dyson index /3 (proven in [9] 
for /3 = 1,2). The centering and scaling parameters are given by 



li w (N,a 1 ,a 2 ) = log 



C+ 



1-6 



o-^(iV,ai,a 2 ) 



iVV(l + + « 2 )(1 + «i + a 2 )C+(l - C+) ' 
Therefore, if Zp ~ F p , then A max ~ 1/(1 + Q-^w-^w^y 

As explained on in [9] (p. 2651) one can write this result informally for A max as 
A max ~ // max + iV- 2/ V max ^ + O (iV" 4 / 3 ) , 

where 

/^max C+ ) 

aLx = [C + (i - C + )] 2 _ 

a/(1 + ai)(l + a 2 )(l + «i + a 2 ) 
Note that this result requires «i,a 2 > 0, so that the case of square or almost square 
matrices is not described by TW. 

By symmetry, one can show that the distribution of A mm (iV, a%, a 2 ) is the same as 
that of 1 — A max (iV, a>2, «i) (note the exchange of ci\ and a 2 ). Therefore, one can easily 
get the typical fluctuations of A m j n from the results presented above. 
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Appendix B. Transformation of Integrals 

The Tricomi solution of the integral equation given by Eq. (11511 can be expanded as 



i _si p K 2 -y/ft-dKCz-*) ^ _«2P f^ 2 v / (*-Ci)(C2-t) ^ 
y;( x ) = 1 2 ^ Jci *(*-*) 2 f_ Jci (l-QCt-x) i- ^ ^ 

n2 V( x -Ci)(C2-x) 

We can further factorize the first integral from (IB.lj) into partial fractions 

fe VC-OKi-t)* = I V(«-C.)(6-«) a _ I f 



i ^(^ — x) X i/^j t — X X J fa t 

Let's start with the first integral in the last expression, where we can use the following 
change of variables: £ = Ci + x — £i + cA, where A = £ 2 — Ci 

Jci t ~ x Jo y~ c 

We can use the following integral 



m'S 9 = W;- C 1, (B.2) 



o J/-c V 2 
where < c < 1, and which leads to 



Similarly, 



Ci 



7T 



where we have used 

1 J til - t) , 7T / , _\ 2 

vv ; dt = -(Vc+T-Vc) . (B.4) 



t + c 2 

Using results (IB .3 1) and ( IB. 41) we can write the above as 
K V(f Cl)(C ^ 
7 Cl i(t - x) x V 

The second integral of (IB.ip can be factorized as 



u (1 — t)(t — x) 1 — Xjfa 1 — t 1 — Xjfa (t — x) 

The second integral in the last expression is the same as before and so we are left with 
the first integral, in which we use the same methods as before to get 



Gathering results, we can finally write the following simplified expression for f{x\ Ci, C2) 

Oil 
X 



27r V(*-Ci)(C2-aO 



(x~ (l-X-y/(l -Cl)(l-C 2 )) +C 
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Since we need a normalized expression for f*(x) i.e. f*(x)dx = 1 (recall that this 
was enforced using the Lagrange multiplier A), this determines C = 2, and thus we have 
the final expression for the limiting spectral density, namely 
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In this appendix we provide details regarding certain integrals which are needed in 
various derivations described in the text. This is done for convenience since some of 
these expressions are not easily available in a single table such as [35]. Some formulas are 
a result of a non-trivial simplification of expressions involving hypergeometric functions, 
and some integrals were performed by the authors. 

Integrals I\ — Xg are used to obtain the simple expression in Equation (|T71) and the 
following shows how this is achieved. After determining the multiplier A 



A = 2ai 1o£ 



—2 + 2a ^{ 2 +21og 4 



where A = ( 2 — Ci -Inserting this into Eq. f JTTj) . we find that S[f*(x); £i> £ 2 ] is given by: 
5[A(x);Ci,Ca] 



-— - y y dxf*(x) logs - — J dxf*(x) log(l 



which can also be expressed as 



cr ai A ai 
S[a,c,A\ = -—- — 



+ a 2 + 2) (log AX! + 22(c)) + «i ^CiG (^X 3 (c) + ^-X 4 ( C ; 



+« 2V / (1-Ci)(l-C 2 ) 0^X5 (a) + ix 6 (a,c) 



A 



A 



do 

-j [(«! + « 2 + 2) (log AXx + X 7 (a)) 



+«iv/C^ f^X 3 (c) + lx 8 (a,cA + « 2V / (1-Ci)(l-C 2 ) f^X 5 (a) + -U 



A 



A 



A 



A 



1-C1 



A' 



and with c > and a > 1. After making the change the variables 



where c 

as in Appendix Appendix B of x = yA + (j, the integrals X& are given by the following 
formulas 



X 2 (c) 
X 3 (c) 
X 4 (c) 
X 5 (a) 
X 6 (a,c) 



v^(l -2/) 
1 log (y + c) 

v^(i - y) 

1 



(y + c) y/y(l - 


y) 


log (y + c) 




{y + c) y/y(i- 


y) 


1 




(a-y) \A/(i - 


y) 


log (y + c) 





dy = Ti 
dy = 27r log 
dy = 
dy = 
dy = 
dy 





+ Vc- 




2 




7T 




C+l) 




71 




C+l) 




71 



log 



4c(c+ 1) (Vc + 1 - Vc 



^a(a — 1) 



(a-y) -z/) 

f log i^faJc+Y) + v^(a - l)c) + log (</a - Va - l) 



ya(a — 1) 
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Mo) 
X 8 (a,c 



o -v) 

log (a - y) 



dy = 2ii log 
dy 



y/a + y/a — 1 



o (y + c) x/y{l-y) 
2tt 



Xq (a) 



v^TT)L 

log (a - y) 



log ( \/a(c+l) + y/(a - l)cj + log (y c + 1 - 



o (a-y) v / ?/( 1 -!/) 



7T 



A/a(a — 1) 



log 



4a(a — 1) (\/a — V& — l) = X 4 (a — 1) 



Gathering these results and after some simple algebra, we eventually get the final 
expression of S[f*(x); Ci, C2] as m equation ([Tj 



(1,(2] = - («! + « 2 + 2) I «i log — ~t — + a 2 log — ^ ^ 2 



+ log v 7 ^ + ?f log - Ci)(l " G) 

+ ai« 2 log log — - — 

To derive the simple expression given in Equation ( l40l) we need the following two 
integrals 

1 

l w (x,a) = 



y/t(l-t)hl(x-t) IT 

t + a 2 



_ y^T) 2 + 2 (2a + 1) In ^ v^+v^^ j 



4 



1—. — , / J (a + 1) x + J a (x - 1) \ 

y Va + 1 + Va y 



(C.l) 



where x > 1 and a > 0, and 
1 

Tu{x,a) = 



y/t(l-t)hi(x-t) dt _n 
t-a 2 







+ Ay/a (a - 1) In 



(^/(a - 1) x + a/ci (x - 1) 
V — 1 + Va 



(v ^_ v^) 2 _ 2 (2a _ 1) , n ^+^T j 

(C.2) 



where x > 1 and a > 1. 

To achieve the simple expression given in Equation (I41j) we need the following two 
integrals 

1 

2i 2 (x,a) = 



y/t(l-t)\n(t + x) 7T 
t+a ' ~ 2 



(v^TT-Vi) 2 + 2(l + 2a)l, 1 ( ^ + 2 vAFTT ) 



— 4 a/ a (a + 



- / y/(a + l)x + y/a (x + l) \ 



(C.3) 
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where x > and a > 0, and 
1 

Xi 3 (x,a) = 



t(l-t)\ll(t + x) _ 7T 
CI L — 

t-a 2 



(y/x + 1 - v^) 2 + 2(1 - 2a) In 



+ 4v ^l) In f ^ Eg± + jT\ 

y Va - 1 + \fa J 



where x > and a > 1. 
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Appendix D. Useful identities for and ( + 

v/(i-C + )(i-C-) 



2 + Oil + «2 
«2 



2 + ai + a 2 

. u u V(l + ai )(l + a 2 )(l + ai + a 2 ) 

^+ ~ ^ ~~ ^Ti i ^ 

(2 + ai + a 2 ) 

^ _ ^^i + 2 + 2a 2 + 2a! + aia 2 

(2 + «i + a 2 ) 2 



/^- /^- 2^1 + «2 

C+ _ V s- 



2 + «i + a 2 ' 

2 v /(l + ai)(l + ai + a 2 ) 
2 + «i + a 2 



2 + Oil + «2 ' 

+ ^ 2v / (1 +„ 2 )( 1 + a . 1 + a2 ) 

2 + ai + a 2 

2 + Cti + a 2 
2sJ\ + a 2 + ai 



VC+(i-C-)->/C-(i-C-4 



2 + CKi + « 2 
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